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ABSTRACT 

Aims. We are carrying out a physical and chemical study of the protostellar envelopes in a representative sample of IM Class 
protostars. In our first paper we determined the physical structure (density-temperature radial profiles) of the protostellar envelopes. 
Here, we study the CO depletion and N2H^ deuteration. 

Methods. We observed the millimeter lines of C^^O, C^^O, N2H+ and N2D+ towards the protostars using the IRAM 30m telescope. 
Based on these observations, we derived the C^^O, N2H^ and N2D^ radial abundance profiles across their envelopes using a radiative 
transfer code. In addition, we modeled the chemistry of the protostellar envelopes. 

Results. All the C^^O 1^0 maps are well fit when assuming that the C^^O abundance decreases inwards within the protostellar 
envelope until the gas and dust reach the CO evaporation temperature, ^^20-25 K, where the CO is released back to the gas phase. 
The N2H^ deuterium fractionation in Class IMs is [N2D^]/[N2H^] =0.005-0.014, two orders of magnitude higher than the elemental 
[D/H] value in the interstellar medium, but a factor of 10 lower than in prestellar clumps. Chemical models account for the C^^O and 
N2H^ observations if we assume the CO abundance is a factor of ~2 lower than the canonical value in the inner envelope. This could 
be the consequence of the CO being converted into CH3OH on the grain surfaces prior to the evaporation and/or the photodissociation 
of CO by the stellar UV radiation. The deuterium fractionation is not fitted by chemical models. This discrepancy is very likely caused 
by the simplicity of our model that assumes spherical geometry and neglects important phenomena like the eff'ect of bipolar outflows 
and UV radiation from the star. More important, the deuterium fractionation is dependent on the ortho-to-para H2 ratio, which is not 
likely to reach the steady- state value in the dynamical time scales of these protostars. 

Key words. Stars: formation - Stars: Emission-Line, Be - Stars: pre-main sequence - Stars: individual: Serpens-FIRS 1 - Stars: 
individual: Cep E-mm - Stars: individual: L1641 S3 MMS 1 - Stars: individual: IC1396N - Stars: individual: CB3 - Stars: individual: 
OMC2-FIR4 - Stars: individual: NGC 7129 FIRS 2 - Stars: individual: S140 - Stars: individual: LkHa 234 



1. Introduction 

Intermediate-mass young stellar objects (IMs) share many char- 
acteristics with high-mass stars (clustering, PDRs) but their 
study presents an important advantage: many are located closer 
to the Sun (d < 1 kpc) and in less complex regions than massive 
star-forming regions. On the other hand, they are also important 
for understanding planet formation since Herbig Ae stars are the 
precursors of Vega-type systems. Despite this, IMs have been 
studied very little so far. A few works on HAEBE stars have 
been carried out at millimeter wavelengths (Fuente et al. 1998, 
2002; Henning et al. 1998), but almost nothing has been done 
for their precursors, the Class IM objects. 

Chemistry has been successfully used to determine the phys- 
ical structure and investigate the formation and evolution of low- 
mass YSOs. Chemical diagnostics have also been shown to be 
good indicators of the protostellar evolution in these objects (see 
e.g. Maret et al. 2004, J0rgensen et al. 2005). However, very few 



works deal with IMs. Fuente et al. (2005a) present a chemical 
study of the envelopes of the Class IM protostar NGC 7129- 
FIRS 2 and the young Herbig Be star LkHof 234. They find that 
the changes in the physical conditions of the envelope during 
its evolution from the Class to the Class I stage (the envelope 
is dispersed and warmed up) strongly influence the molecular 
chemistry. The Class object NGC 7129-FIRS 2 presented evi- 
dence of H^^CO^ depletion. Moreover, the deuterium fractiona- 
tion, measured as the DCO^/H^^CO^ ratio, decreases by a factor 
of 4 from the Class to the Herbig Be star, very likely owing 
to the increase in the kinetic temperature. Regarding the abun- 
dance of complex molecules, the beam-averaged abundances of 
CH3OH and H2CO increase from the Class to the Herbig Be 
star. A hot core was also detected in NGC 7129-FIRS 2 (Fuente 
et al. 2005 a,b). Although two objects are not enough to establish 
firm conclusions, these pioneering results suggest that chemistry 
is also a good indicator of the evolution of IMs. 
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We are carrying out a chemical study of a representative sam- 
ple of IM Class YSOs This is the first systematic chemical 
study of IM Class objects that has been carried out so far. 
Some properties like the temperature of the protostellar enve- 
lope and the clustering degree depend on the final stellar mass, 
so the results for low-mass stars cannot be directly extrapo- 
lated to intermediate-mass objects. In the first paper (Crimier 
et al. 2010, hereafter CIO), we determined the physical struc- 
ture (density-temperature radial profiles) by modeling the dust 
continuum emission. We now present the observations of the 
millimeter lines of C^^O, C^^O, N2H+, and N2D+ in the same 
sample. Our goal is to investigate the CO depletion and N2H^ 
deuteration in these Class YSOs. For comparison, we also in- 
clude 2 Class I objects, LkHa 234, and S140. 



2. Observational strategy 

Our selection was made to have a representative sample of Class 
IM YSOs, including targets with diff'erent luminosities (40 - 
10^ Lq) and evolutionary stages. An important complication in 
the study of massive stars is that they are located in complex 
regions and are therefore difficult to model. The targets in this 
sample were chosen to lie preferentially in isolated areas with 
respect to the 30m telescope beam. We also selected sources for 
which continuum maps at submillimeter and/or millimeter wave- 
lengths are available in order to be able to model their envelopes. 
The list of sources and their coordinates are shown in Table 1. 
To provide a comparison with Class I sources, we added SI 40 
and LkHof 234 to the sample. 

Most of the observations reported here were carried out with 
the IRAM 30m telescope at Pico de Veleta (Spain) during three 
diff'erent observing periods in June 2004 in position switching 
mode. Our strategy was to first make long integration single- 
pointing observations towards the star position and then to make 
96'' X 96" maps around the center position. The maps were sam- 
pled with a spacing of 12'' in the inner 48" regions and 24" 
outside. The only exceptions were LI 641 S3 MMSl and S140. 
In L1641 S3 MMSl, we only observed a radial strip in C^^O 
and C^^O. We did not observe C^^O and C^^O maps in S140. A 
summary of the observations for each source is shown in Table 
1, and the list of observed lines and the telescope characteris- 
tics is shown in Table 2. During the observations, lines of the 
same species were observed simultaneously using the multire- 
ceiver capability of the 30m telescope. In this way, we mini- 
mized relative pointing and calibration errors. As backends we 
used in parallel an autocorrelator split into several parts provid- 
ing a spectral resolution that was always better than ~78 kHz and 
a 1 MHz-channel-filter-bank. Examples of the single-pointing 
observations are shown in Figs. 1 and 2. The intensity scale is 
the main brightness temperature. 

Observations of the N2H+ J=4^3 line (freq=372.6725 GHz) 
were carried out towards Cep E-mm, IC 1396 N, NGC 7129- 
FIRS 2, and L1641 S3 MMSl using the JCMT telescope at the 
Mauna Kea (Hawaii). In all these sources, we carried out small 
maps of 75"x75" with a spacing of 25", but the emission was 
only detected towards the star position (see Fig. 3). All the lines 
were observed with a spectral resolution of 0.488 MHz. The 
spectra of the J=4^3 line are shown in Fig. 3. 

In this paper, we also use the C^^O, N2H+ and N2D+ data 
towards NGC 7129-FIRS 2 and LkUa 234, which has already 
been published by Fuente et al. (2005b). 



Table 2. Description of the observations 



Line 




HPRW 
nrD vv 




ICl 


C^^O 1^0 


109782.2 


22" 


0.75 


IRAM 




219560.3 


11" 


0.55 


IRAM 


C^^O 1^0 


112359.3 


22" 


0.74 


IRAM 




224714.4 


11" 


0.54 


IRAM 


N2H+ 1^0 


93173.2 


26" 


0.77 


IRAM 


N2H+ 4^3 


372672.5 


13" 


0.63 


JCMT 


N2D+ 2^1 


154217.1 


16" 


0.67 


IRAM 


N2D+ 3^2 


231321.7 


11" 


0.52 


IRAM 



3. Results 

The spectra of the C^^O 1^0, C^^O 1^0, N2H+ 1^0 and N2D+ 
2^1 lines towards the star position are shown in Figs. 1 and 2. 
The integrated intensity maps of the C^^O 1^0, N2H^ 1^0, and 
N2D^ 2^1 lines are shown in Figs. 4 (Class 0) and 5 (Class I). 
In all the Class sources, the emission of the N2H^ 1 ^0 line is 
compact and peaks towards the star position, revealing a dense 
core around the star. However, towards the Class I sources, the 
peaks of the N2H^ emission are offset from the far-IR source. 
For the Herbig Be star LkHa 234, the emission peak of the N2H^ 
1^0 line is located at an offset (-6", 18") from the star position 
(see Fuente et al. 2005a and Fig. 5). In S140, the N2H+ emission 
peaks in an arc-shaped feature surrounding the sources IRS 1, 
2 and 3. In Class I sources, either the N2H^ is destroyed by the 
evaporated CO or the outflow, or the UV radiation from the star 
has already disrupted the parent core. 

The emission in the C^^O 1^0 line usually presents a dif- 
ferent morphology from that of N2H^ 1^0 line. In CB3 and 
IC 1396N, the emission in the C^^O line has an elongated 
shape, much more extended than that of N2H"^ (see Fig. 4). 
In 0MC2 FIR 4 and NGC 7129-FIRS 2, the emission of the 
C^^O line surrounds the star position instead of having a max- 
imum towards it (see also Fuente et al. 2005a). Only in the 
low-luminosity sources Serpens-FIRS 1 and Cep E-mm does the 
emission from the C^^O and C^^O lines peak at the star position 
and present a morphology similar to that of the N2H^ emission. 

For the sources in which the signal-to-noise ratio of the 
N2D+ map is high enough, Serpens-FIRS 1 and NGC 7129 
FIRS 2, we compared the N2D+ 2^1 and N2H+ 1^0 maps. In 
Serpens-FIRS 1, the N2D"^ 2^1 emission does not peak towards 
the star, but does in the case of NGC7129-FIRS 2. 

In Fig. 6, we show the mean integrated intensity emission 
of the C^^O 1^0, N2H+ 1^0, and N2D+ 2^1 lines in con- 
centric rings around the protostar. In all the sources, the emis- 
sion of the N2H^ 1^0 line decreases outwards from the star. 
The radial profiles of the C^^O 1^0 line, however, greatly differ 
from one source to the next. The presumably youngest sources, 
OMC2 FIR 4 and NGC 7129-FIRS 2, show a flat profile. This is 
the expected picture when the CO abundance decreases towards 
the center, mainly because the molecules are frozen onto the 
grain surfaces. In Serpens-FIRS 1, Cep E-mm, IC 1396 N, and 
CB3, the C^^O 1^0 emission decrease with the distance from 
the star, similar to N2H"^ but with a less steep profile. Intense 
emission of the C^^O 1^0 line is detected at the edges of the 
protostar envelopes, which shows that a lower density envelope 
also contributes to the emission of this line. This envelope con- 
tribution is especially important for Serpens-FIRS 1 . In fact, we 
must model the envelope contribution in order to fit the C^^O 
1^0 emission from this protostellar core (see Appendix Al). 
One could interpret the different radial profiles of the C^^O 1^0 



T. Alonso-Albi et al.: Chemical study of intermediate-mass (IM) Class protostars. 

Table 1. Selected sample 



3 



Object 


RA(2000) 


Dec(2000) 


Lum. (Lq) 


d 




Serpens-FIRS 1 


18:29:49.6 


+01:15:20.6 


33 


230 


Maps in C^^0,C^^0,N2H+ and N2D+ 


Cep E-mm 


23:03:13.1 


+61:42:26.0 


100 


730 


Maps in Ci^O,Ci80,N2H+ and N2D+ 


L1641 S3 MMS 1 


05:39:55.9 


-07:30:28.0 


67 


500 


Strip in C^^O and C^^O, maps in N2H+ and N2D+ 


IC1396N 


21:40:41.7 


+58:16:12.8 


150 


750 


Maps in Ci^0,Ci^0,N2H+ and N2D+ 


CB3 


00:28:42.7 


+56:42:07.1 


1000 


2500 


Maps in C^^CCi^CNsH^ and N2D+ 


OMC2-FIR4 


05:35:26.7 


-05:10:00.5 


1000 


450 


Maps in Ci^O,Ci80,N2H+ and N2D+ 


NGC 7129-FIRS 2 


21:43:01.7 


+66:03:23.6 


500 


1250 


Strip in C^^O, maps in C^^0,N2H+ and N2D+ 


S140 


22:19:18. 1 


+63:18:54.6 


10^ 


910 


N2H+ and N2D+ maps 


LkHor 234 


21:43:06.8 


+66:06:54.4 


500 


1250 


Maps in Ci^O,N2H+ and N2D+ 



Tmb(K) 



Serpens — FIRS1 



Cep E— mm 



L1641 S3 mm 



IC 1396 N 



: c'^o 1-0 






_i 1 1 1 1 1 1 1 1 1 
: c^^o 2-1 


1 1 1 1 1 1 1 1 i_ 
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1 NgH"^ 1-0 h 
:l 1 mTTTTi I 




~ NjO"" 2-1 
^^^^^^ 
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5 
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-5 5 
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Fig. 1. Examples of the single-pointing observations towards the Serpens-FIRS 1, Cep E-mm, LI 641 S3 MMSl, and IC 1396 N. 
All the spectra were observed with the 30m telescope. The intensity scale is the main brightness temperature. 
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Fig. 2. The same as Fig. 1 for CBS, 0MC2 FIR 4, NGC 7129-FIRS 2, and S140. 
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^1 11 1 1 , 11 1 1 

I NGC 7129 


N2H + (4-3) ; 
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-15 -10 -5 
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3 
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-5 

V|sr(km 


5 
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Fig. 3. Spectra of the N2H+ 4^3 line towards NGC 7129- 
FIRS 2, Cep E-mm, IC 1396 N, and L1641 S3 MMSl. 



emission as an evolutionary trend, with the youngest protostars 
having flat profiles while the oldest have steep profiles. Given 
the complexity of these regions, however, one should be cau- 
tious and use multiple evolutionary tracers to define relative age. 
For instance, a lack of C^^O emission towards the star position 
could come from CO depletion in the case of a very young object 
or to photodissociation for a borderline Class O/I. The radial pro- 
files of the N2D+ 2^1 (3^2 for NGC 7129-FIRS 2) are more 
uncertain because of the low S/N ratio of the maps. 



4. LTE column densities 

We have derived the C^^O column densities using the rotation 
diagram method. This method gives an accurate estimate of the 
column density provided that the emission is optically thin, is 
thermalized, and arises from a homogeneous and isothermal 
slab. In the case of a density and temperature distribution, the 
derived column density represents an average value over the ob- 
servational beam. 

The excitation of C^^O and C^^O are very similar, and we 
have better signal-to-noise spectra for the more intense C^^O 
lines. For this reason, we derived the C^^O column density by 
assuming optically thin emission and the rotation temperature 
derived from the C^^O data. The molecule N2H^ presents hy- 
perfine splitting. This allows us to estimate the line opacity di- 
rectly from the hyperfine line ratios. We derived the total N2H^ 
column density from the opacity of the N2H^ 1^0 line, and as- 
sumed the rotation temperatures derived from N2D^ data when 
there was no N2H^ 4^3 observations. The same rotation tem- 
perature was always used for both N2H^ and N2D^. We assumed 
a beam filling factor of 1 for all the lines regardless of the ob- 
servational beam size. This approach overestimates the rotation 
temperature when the emission is centrally peaked. In Table 3 we 
present the derived C^^O, C^^O, N2H+, and N2D+ column densi- 
ties. We also show the NCC^^OVNCC^^O) (hereafter Ri) and the 
N(N2D+)/N(N2H+) (hereafter R2) ratios. For the whole sample, 
we obtain a value for Ri around 3.3+/-0.3, the expected value for 
optically thin emission. In the worst case, Ri=2.5, the correction 
to the C^^O column density because of the opacity is only a fac- 
tor -1.3. Thus, opacity eff'ects are not important in our C^^O 
column density estimates. 



4.1. /V2D+ andN2H^ 

The deuterium fractionation of N2H^ (R2) strongly depends on 
the CO depletion factor and the gas temperature (Caselli et al. 
2002, CeccarelH & Dominik 2005, Daniel et al. 2007). In our 
sample of IM Class protostars, we derive values of R2 rang- 
ing from 0.005 to 0.014. These values are 3 orders of magnitude 
higher than 10"^, the elemental value in the interstellar medium 
(Oliveira et al. 2003). Nevertheless, the R2 values are a factor 
of 10 lower than those found in prestellar clumps by Crapsi et 
al. (2005). According to the values of R2, we can classify our 
sources into two groups: (i) highly deuterated sources that have 
values of R2>0.01 and (ii) moderately deuterated sources with 
R2<0.01. Cep E-mm, CB3, and NGC 7129-FIRS 2 belong to 
the first group. Assuming that the N(N2D^)/N(N2H^) ratio is a 
good gas temperature indicator, these protostars should be the 
coldest and very likely the youngest of our sample. We have to be 
cautious with CB3, however. Since this source is the most distant 
(d=2500 pc), our single-pointing observations trace a larger frac- 
tion of the envelope. The sources Serpens-FIRS 1, IC 1396 N, 
0MC2-FIR 4, and S140 belong to the second group. S140 is 
a Class I YSO and IC 1396 N is considered a borderline Class 
O/I YSO. The results are thus consistent with our interpretation 
of the objects in this group as being more evolved. In Serpens- 
FIRS 1, however, the N2D^ 1^0 emission comes mainly from 
the lower density envelope and the ratio cannot be considered a 
tracer of the evolutionary stage of the protostar. 0MC2 FIR 4 
is also difficult to classify. The deuterium fractionation suggests 
an evolved object but the morphology of the C^^O map is more 
consistent with a young Class star. The peculiarity of this pro- 
tostellar envelope has already been pointed out by Grimier et al. 
(2009). They find that the envelope of 0MC2 FIR 4 is pecuHarly 
flat and warm with a radial density power-law index of 0.6. 

We present the N2H^ deuterium fractionation as a func- 
tion of the N(C^^0)/N(N2H+) ratio in Fig. 7. An excellent 
correlation between these two quantities is found in prestellar 
cores with the deuterium fractionation decreasing with increas- 
ing N(C^^0)/N(N2H+) ratio (see Crapsi et al. 2005). This cor- 
relation is not valid for IM Class protostars. As discussed in 
the following sections, this is due to the complexity of these 
intermediate mass star-forming regions. In prestellar cores the 
chemistry of these species is only driven by the CO depletion. 
In IMs, other phenomena like photodissociation and shocks are 
also playing important roles. 

5. Radiative transfer model 

We utilized a general ray-tracing radiative transfer cod^Jo de- 
rive the fractional abundance profiles of C^^O, N2H^ and N2D^ 
across the envelopes. Assuming appropriate radial profiles for 
the temperature, density, molecular abundance, and turbulence 
velocity, this model calculated the brightness temperature dis- 
tribution on the sky. The model map was then convolved with 
the telescope beam profile. The underlying source geometry was 
assumed to be a sphere, with the inner and outer radii, and 
temperature-density (T-n) profiles derived by CIO. The size of 
the grid was set to 32x32 cells. The cells have diff'erent sizes 
along the line of sight to account for the diff'erent slopes in the 
density and temperature profiles. We used very small cells (<100 
AU) near the center, where the temperature and density gradients 
are highest. In the outer regions, the cells reach several thousand 



^ The model is called DataCube and a link to install it is available 
upon request. 
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Continuum 850 /xm C^^O 1-0 N2H^ 1-0 N2D^ 2-1 




40 20 -20 -40 

Aa (arcseconds) 

■ # T ▲ • 

Main 24 jam Secondary 24 /xm 450 /xm 850 Peak N2H 30m 



Fig. 4. Continuum maps at 850 yum from CIO (left column), and integrated intensity maps of the C^^O 1^0, N2H^ 1^0 and N2D^ 
2^1 (3^2 in the case of NGC 7129-FIRS 2) lines towards the observed Class IMs. In each panel, we have drawn the area 
and the number of rings considered in our model fitting. The emission peak at several wavelengths is shown with diff'erent marks. 
Contour levels starts at 3cr level and are a. 0.5, 1.6, 2.7, 3.8, 4.9, and 5.4 Jy/beam; b. 1 to 9 K km s"^ by steps of 1 K km s"^ ; c. 3 to 
9 K km s-i by steps of 3 K km s'^ d. 0.5 to 1.5 K km s'^ by steps of 0.5 K km s'^ e. 0.16, 0.49, 0.81, 1.14, and 1.46 Jy/beam; f. 

0. 5 to 3 K km s"^ by steps of 0.5 K km s'^ g. 0.5 to 2.5 K km s'^ by steps of 0.5 K km s'^ h. 0.3, 0.9, 1.5, 2.1, and 2.8 Jy/beam; 

1. 2.0 to 8.0 K km s'^ by steps of 2.0 K km s'^ j. 3.0 to 9.0 K km s'^ by steps of 3.0 K km s'^ k. 0.07, 0.20, 0.33, 0.46, and 0.59 
Jy/beam; 1. 2.0, 4.0 K km s'^ m. 1.0 to 3.0 K km s"^ by steps of 1.0 K km s'^ n. 1.0, 2.5, 4.0, 5.5, and 7.0 Jy/beam; o. 0.5 to 
3.5 K km s"^ by steps of 0.5 K km s'^ p. 4.0 to 24.0 K km s"^ by steps of 4.0 K km s'^ q. 2.1, 6.3, 10.6, 14.8, and 19.0 Jy/beam; 
r. 2.0 K km s'^ ; s. 3.0, 6.0 K km s'^ ; t. 0.5 K km s'^ . 
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LkH« 234 



S140 



40 20 -20 -40 

Aa (arcseconds) 



Fig. 5. Integrated intensity maps of the C^^O 1^0 and N2H^ 1^0 lines towards LkHa 234 and S140. Contour levels are: a. 0.07, 
0.15, 0.23, 0.31, 0.39, 0.46, 0.54, 0.62, and 0.69 Jy/beam; b. 2 to 5 K km s"^ by steps of 1 K km s'^ c. 1 to 5 K km s"^ in steps of 
1 K km s"^ ; d. 1 to 8.0 Jy/beam in steps of 1 Jy/beam ; e. 3 to 21 K km s"^ in steps of 3 K km s"^ 

Table 3. LTE column densities 



Source 


T. 


N(C^^O) 


N(C^^O) 


RJ 


T, (K) 


N(N2H+) 


N(N2D+) 


R* 




(K) 


(cm-2) 


(cm-2) 




(K) 


(cm-2) 


(cm-2) 




Serpens-FIRS 1 


11 


8.2x10^^ 


3.2x10^^ 


2.6 


11^ 


7.7x101-^ 


3.4x1011 


0.004 


Cep E-mm 


15 


2.5x10^5 


7.6x1014 


3.3 


gb 


1.4x1013 


5.8x1011 


0.04 


L1641 S3 MMS 1 


10 


3.2x10^5 


1.3x10^5 


2.5 


8^ 


3.6x1013 


<2.1xl0ii 


<0.006 


IC1396N 


25 


1.1x10^6 


3.7x10^5 


3.0 


10^ 


41x1013 


<41xl0ii 


<0.01 


CB3 


14 


5.2x10^5 


1.6x10^5 


3.3 


6^ 


1.2x1013 


1.0x1012 


0.08 


OMC2-FIR4 


22 


48x10^5 


1.2x10^5 


34 


22^ 


2.0x1014 


8.7x1011 


0.004 


NGC7129-FIRS2 


23 


1.9x10^5 


6.9x1014 


2.8 


13" 


3.8x1013 


5.4x1011 


0.014 


S140 


50^ 


1.0x10^6 


2.7x1015 


3.7 


50^ 


1.2x1014 


<4.5xl0ii 


<0.004 


LkHof 234 ^ 










22 


1.0x1013 


<2.3xl0ii 


<0.02 



* Ri=N(Ci^O)/N(Ci^O), R2=N(N2D+)/N(N2H+). 

^ Rotation temperature derived from our N2D^ observations. 

^ Rotation temperature derived from our N2H+ observations. 

' From Fuente et al. (2005a). 

^ From Minchin et al. (1995). 



AU in size. The turbulent velocity was assumed to be fixed at 
1.5 km s"\ consistent with the linewidths of the observed lines. 
Finally, in each cell, the excitation temperature was calculated 
with the RADEX code (van der Tak et al. 2007), which uses the 
slab LVG approximation at each shell. We used the collisional 
rates provided by the LAMDA databasfl (Schoier et al. 2005). 
For N2D^, we used the same collisional rates as for N2H^ . 



^ LAMDA database is available at 

|http ://w w w. st rw.leidenuniv.nl/~moldata/[ 



With these assumptions, we searched for the best fit on every 
source and molecule observed. The fitting process consisted of 
averaging the observed and modeled fluxes in concentric rings 
around the center position using GILDAS software. The cen- 
ter position was selected to be the continuum emission peak at 
850 yum. The radius of the first ring was set to HPBW/4, and in- 
cremented by the same amount in consecutive rings. As a first 
step we searched for the best fit using a constant abundance. 
Since this approach seldom produced good fits, we next searched 
for the best fit using radial power functions and step functions. 
These functions were not selected arbitrarily but were selected 
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NTENSITY VARIATION WITH RADIUS 




I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I 
10 20 30 10 20 30 10 20 30 




10 20 30 10 20 30 1 20 30 

• C^'O (1-0) r ["] 

♦ N2H^ (1-0) 

■ N2D^ (2-1) [(3-2) in N7129] 

Fig. 6. Radially integrated intensity profiles of the C^^O 1^0, N2H+ 1^0, and N2D+ 2^1 (3^2 for NGC 7129-FIRS 2) lines in 
our sample of Class IMs. Each point represents the mean value of the line integrated emission in a ring of a given radius. The 
errors are the rms of the line integrated emission in the ring. The step between two consecutive points, which is equal to the size of 
each ring, is equivalent to 1/4 of the beam at each frequency. 



Table 4. Temperature-density (T-n) profiles from Grimier et al. 2010 



Source 


CBS 


Cep-E 


IC 1396 N 


NGC 7129 


Serpens 




mm 


mm 


B1MA2 


FIRS 2 


FIRS 1 


RA(J2000) 


00:28:42.1 


23:03:12.7 


21:40:41.8 


21:43:01.5 


18:29:49.8 


Dec(J2000) 


56:41:59.4 


61:42:27.4 


58:16:13.5 


66:03:25.0 


01:15:18.4 


Dust optical depth at 100//m, tioo 


5.8 


5.0 


1.4 


2.3 


3.0 


Density power law index, a 


2.2 


1.9 


1.2 


1.4 


1.5 


Envelope thickness, ToutM 


400 


500 


630 


180 


200 


Inner envelope radius, rin, (AU) 


260 


70 


50 


100 


30 


Outer envelope radius, rout, (AU) 


103000 


35800 


29600 


18600 


5900 


Radius at T^ust = 100 K, tiook, (AU) 


700 


223 


180 


373 


102 


H2 density at riooK, ^0, (cm~^) 


7.5x10^ 


2.0x10^ 


4.3x10^ 


4.4x10^ 


2.2x10^ 


Envelope mass, Mgnv, (Mq) 


120 


35 


90 


50 


5.0 



to mimic the predictions of chemical models. The step func- 
tion accounts for the abrupt sublimation of the GO ices thanks 
to the increase in the dust temperature going inwards, whereas 
the power-law profile accounts for a smoother change in the 
abundance of the species. The angular resolution of our observa- 
tions (~16"-27" depending on the frequency, i.e, -16000 AU- 
27000 AU at a distance of 1000 pc) prevents us from tracing the 
inner region (R< a few 1000 AU) of the protostellar envelope. 
For this reason, we assume a constant abundance in the inner 
part. 

We fit the integrated intensity maps of the G^^O 1^0, N2H+ 
1^0, N2H+ 4^3, and N2D+ 2^1 lines, when possible. The 
maps of the G^^O 2^1 and G^^O 1^0 and 2^1 lines are of less 
quality so we did not use them in our fits. Some sources show 
deviations from the spherical symmetry because of the contribu- 
tion of the surrounding molecular cloud (Serpens-FIRS 1, GB3) 



or because of the cavity excavated by the outflow (IG 1396 N). 
In these cases we masked part of the map to avoid this contam- 
ination in the rings. A short description of earch source and the 
individual details about the modeling are given in Appendix A. 
The areas masked are shown in Fig. 4, and the abundance pro- 
files obtained from the model are listed in Table 5. We did not 
model the Glass I sources S140 and LkHa 234. The modeling of 
0MG2 FIR 4 will be the subject of a forthcoming paper. 



6. Limitations of our model: CO evaporation 
temperature 

The model used here is not a reliable predictor of the chemical 
and physical properties within the inner envelope for several rea- 
sons. First of all we only considered the low- J rotational molec- 
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Table 5. Abundance profiles derived from the observations 



Serpens-FIRS 1 


Xo 


Ro (AU) 


Xi 


Ri (AU) 


X2 


RMS (K km s" 


1) 


Diameter(")/HPBW(") 


C^^O 
N2H+ 


1.4x10"' 
4.0x10"!^ 


2000 
3000 


<7xl0"^ 
<3xlO"ii 


6000 
6000 


1.1x10"' 
2.5x10-10 


0.27 
0.29 




2.3 
1.9 


N2D+ 




<2xl0-i^ 




6000 


1.9x10-11 


0.6 




3.1 


t:, (k) 

Kv (cm ^) 


20 

2.4x10^ 
















Cep E-mm 


Xo 


Ro (AU) 


Xi 


Ri (AU) 


X2 


RMS (K km s" 


1) 




C^^O 

N2H+ 


6.0x10-^ 
2.4x10-10 


3500 
6000 


l.lxlO-'x(r/35800) 
2.5xlO-iOx(r/35800) 


0.24 
0.28 




4.4 
3.7 


N2D+ 






7.0x10"!^ 










6.1 


T.V (K) 

Yley (cm-3) 


20 

1.0x10^ 
















IC 1396 N 


Xo 


Ro (AU) 


Xi 


Ri (AU) 


X2 


RMS (Kkm^- 






C^^O 


6.0x10-^ 


5500 


1.25xl0-'x(r/29600) 


0.24 




3.5 


N2H+ 


4.2xl0"i" 


10000 


1.0xlO"ii 


15000 


1.8xlO"io 


0.50 




3.0 


N2D+ 






<7xl0"i^ 










4.9 


T.V (K) 
(cm ^) 


19 

7.0x10^ 
















CB3 


Xo 


Ro (AU) 


Xi 


Ri (AU) 


X2 


RMS (K km s" 


1) 






1.3x10"' 


25000 


<3.0xlO"« 


60000 


9.0x10"' 


0.24 




3.7 


N2H+ 
N2D+ 






6.5xlO"i" 
3.0xlO"ii 






0.45 




3.1 
5.1 


Tev (K) 

(cm-3) 


15 

3.0x10^ 
















NGC 7129-FIRS 2 


Xo 


Ro (AU) 


Xi 


Ri (AU) 


X2 


RMS (K km s" 


') 




C^^O 
N2H+ 






4.0x10-^ 
4.7xlO"io 






0.4 
0.4 




1.4 
1.1 


N2D+ 


<3.0xlO-i^ 


9000 


4.5xl0"ii 


11000 


<3.0xlO"i^ 


0.13 




2.7 



* Temperature and density at the CO evaporation radius. 



ular lines. The emission of these transitions arises mainly from 
the outer envelope. This fact, together with the limited spatial 
resolution of our single dish observations, prevent us from de- 
termining the variation in the abundance in the inner envelope. 
Additionally, in our modeling we use the n-T profiles derived 
by CIO. These profiles are a reasonable approximation for the 
outer envelope but are relatively unconstrained in the inner re- 
gion. Some parameters, like the inner radius of the envelope and 
the dust temperature at this radius, are thus poorly determined. 

Another important source of uncertainty is the degeneracy 
of the solutions. For instance, in the case of C^^O we can obtain 
similar results by varying the molecular abundance in the inner 
region (Xq) or the CO evaporation radius (Ro), as long as the 
line is optically thin and the number of molecules in the beam 
remains constant. Here we discuss the impact of this degeneracy 
in the derived CO evaporation radius. 

We ran a grid of models for Cep E-mm, IC 1396 N, and 
CBS. These are the sources in which the ratio between the angu- 
lar diameter of the envelope and the HPBW of the telescope is 
greatest, so we are more sensitive to the spatial variations of the 
chemical abundances. In all these models, we assume that the 
C^^O abundance profile is X=Xo for radii less than the evapo- 
ration radius, Rq, and X= XoX(R/RoutT for radii larger than Rq. 
For each value of Xq, we fit Rq and a. Xq is the un-depleted 
C^^O abundance, Rq the evaporation radius, and a measures the 
gradient in the C^^O abundance due to depletion. The canonical 
abundance of C^^O in principle depends on the Galactocentric 
distance (DGC). Following Wilson & Matteucci (1992), the CO 
abundance is given by 

X^co) = 9.5 x 10^exp(lA05 - (0.13 x DGC(kpc))). (1) 



Following Wilson & Rood (1994), the oxygen isotope ra- 
tio ^^O/^^O depends on DGC according to the relationship 
160/I80 = 58.8 x DGC(kpc) -h 37.1. At the distance of the 
Sun, 8.5 kpc, Xci8o= 1-7x10"^. This is an average value, so 
there may be local effects. In our models we varied Xq from 
1x10"^ to 4x10"^, the range of Xq values that we consider 
reasonable. For each value of Xq, we varied Rq and a to find 
the best fit. In Fig. 8 we show the results of our models for 
Cep E. The best fit is Xq = 1x10"^ Ro=2510 AU, a=l with 
RMS=0.37 K km s"^ We consider that the models with RMS 
values departing by a factor of less than 1.3 from the mini- 
mum value are still acceptable. Following this criterion, we have 
two different families of solutions. Assuming Xo=lxlO"^, the 
best fit is obtained with a^l and Ro~2100-3100 AU. Assuming 
Xo=2xlO"^, we have reasonably good solutions with and 
Ro -2200-3200 AU. We have not found any acceptable solution 
for higher values of Xq. The evaporation radius of C^^O is thus 
2600+600 AU. These radii correspond to dust temperatures of 
20-25 K. Interestingly, we would have a better fit to the obser- 
vations (RMS=0.24 K km s"^) if we allowed the C^^O abun- 
dance in the inner region to fall below 1x10"^. This is the solu- 
tion we show in Table 5. We do not reject this solution because 
there are several physical reasons that the CO abundance could 
be lower in the inner envelope. First of all, C^^O could have 
been photodissociated in the regions very close to the recently 
formed star. Alternatively, part of the CO could have been trans- 
formed into CH3OH on the icy mantles, and the CO abundance 
once evaporated would then be different from the initial value. 
However, given the large uncertainty in the physical structure of 
the inner parts of the envelope, we cannot draw any conclusion. 
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In Fig. 9 we show the results for IC 1396 N. The 
best fit is Xo = 1x10"^ Ro=2000 AU, and a=0.5 with 
RMS=0.36 K km s~^ . We consider that the models with values of 
RMS<0.5 are acceptable. Following this criterion, we only find 
solutions for Xo=lxlO"^ and Ro<3000 AU; i.e., the CO evapo- 
ration temperature should be >23 K. Again, the fit can be sig- 
nificantly improved if values of Xq lower than 1x10"^ are con- 
sidered. This is similar to the case of Cep E. Finally, in Fig. 10, 
we show the same grid of models for CB3. The best fit is for Xq 
= 4x10"^ Ro=4000 AU, and a=0J6 with RMS=2.0 K km s-\ 
which corresponds to a CO evaporation temperature of ~30 K. 
However, this is not a good fit. In fact, we obtain a much better 
solution when assuming a step function and Tev~15 K (see Table 
5). 

For Serpens-FIRS 1 and NGC 7129 FIRS2, the protostars 
are barely resolved by our single-dish observations. In addition, 
the contribution of the surrounding cloud is very significant. For 
these reasons, we did not carry out a study similar to what is 
described above. We looked only for abundance profiles that fit 
our observations and are consistent with the behavior predicted 
by chemical models (see Table 5). In these two cases, we can 
fit the observations when assuming that the CO is evaporated at 
temperatures around 20-25 K. We conclude, therefore, that our 
C^^O observations towards the Class stars are better fit assum- 
ing CO evaporation temperatures of 20-25 K that are consistent 
with the CO being bound in a CO-CO matrix. 

The fits to the N2H^ and N2D^ emission profiles have been 
done by hand. For N2H^ we tried two options :(i) constant N2H^ 
abundance and (ii) an abundance profile with a radial variation 
similar to that of C^^O. Option (ii) follows the expectation that 
the N2H"^ abundance is strongly dependent on the CO abun- 
dance. In the case of CB3 and NGC 7129 FIRS 2, we are able 
to fit the N2H^ observations by assuming a constant abundance. 
For the rest of sources, the N2H^ observations are better fit by 
assuming that the N2H^ is depleted in the cold regions similarly 
to CO. 

One expects an annular distribution for the N2D^ abundance, 
with the maximum N2D^ abundance in the region with the great- 
est CO depletion. The morphology observed in the integrated 
intensity map of the N2D^ 2^1 line towards Serpens-FIRS 1 
suggests that the N2D^ emission is dominated by the foreground 
molecular cloud. For this reason we selected a one-step func- 
tion with a constant abundance inside and outside the protostel- 
lar core. In the case of NGC 7129-FIRS 2, we detected a clump 
of N2D^ 3^2 emission towards the star. However, interferomet- 
ric observations show that on scales of a few 1000 AU, there is a 
lack of emission towards the star (Fuente et al 2005a). We thus fit 
the N2D^ emission assuming an annular abundance distribution. 

7. Chemical model 

The chemical composition was modeled with the simple chem- 
ical code originally described in Caselli et al. (2002), and up- 
dated following Caselli et al. (2008) with new measurements of 
the CO and N2 binding energies (Collings et al. 2003; Oberg et 
al. 2005) and sticking coefficients (Bisschop et al. 2006), ther- 
mal desorption, and the detailed physical structure. The clouds 
are assumed to be spherically symmetric, with the density and 
temperature profiles derived from the dust continuum emission. 
An interpolation procedure has been included in the code in 
order to have smooth profiles. The chemical network contains 
the neutral species CO and N2, which can freeze out onto dust 
grains and desorb owing to cosmic ray impulsive heating (as 
in Hasegawa & Herbst 1993) and by thermal evaporation (fol- 



NGC 7129 FIRS 2 



0MC2 FIR 4 



Li 641 S3 MMSl 



S140 V" Serpens FIRS 1 



N(C^^O)/N(NgH"') 



Fig. 7. Deuterium fractionation of N2H^(R2) vs the 
N(C^^0)/N(N2H+) ratio in the studied YSOs. Inverted empty 
triangles are upper limits for L1641 S3 MMSl, IC 1396 N, and 
S140, where we did not detect the N2D+ 2^1 line. 



CepE error surface (RMS (K km/s)) CepE 



face (RMS (K km/s)) 




1000 2000 
Ro (AU) 



1000 2000 
Ro (AU) 



Fig. 8. Plots of the RMS, defined as ^{Imodd - hbsf) for the grid 
of models run to reproduce the C^^O emission in Cep E-mm. 
The C^^O abundance is assumed to be X=Xo for radii less than 
Ro, and X= XoX(R/Roi^O^ for radii larger than Rq. The values of 
Xo are 1.0x10"^ (panel a.), 2.0x10"^ (panel b.), 3.0x10"^ (panel 
c), and 4.0x10"^ (panel d.). The solution with the lowest RMS, 
0.37, is indicated by a cross. We consider acceptable solutions 
those with RMS less than 0.5. Contour levels are 0.5, 1, 2, and 5 
Kkms-^ 



lowing Hasegawa et al. 1992). The initial abundances of CO 
and N2 have been fixed to 9.5x10"^ (Frerking et al. 1982) and 
2x10"^, respectively. The abundances of molecular and atomic 
nitrogen are difficult to determine in dense cores, but recent 
works suggest low values in low-mass, star-forming regions: 
X(N) = n(N)/n(H2) < 2x10"^ (Hily-Blant et al. 2010), and 
X(N2)~10"^ (Maret et al. 2006). Our adopted value is consistent 
with 25% of the total abundance of N (n(N)/n(H) = 7.9x10"^ 
see Anders & Grevesse 1989) locked in N2 and with the pseudo- 
time dependent models of Lee et al. (1996), after the chem- 
istry reaches steady state. Although atomic oxygen can aff'ect 
the amount of deuterium fractionation (see discussion in Caselli 
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Table 6. Model results 



Source 


Model 


CO Binding energy 




XN2H+ 


XN2D+ 




Serpens-FIRS 1 


Model 3 


1100 


431 


3.89 


4.52 


only 10% of CO evaporated 




Model 4 


1100 


2.55 


3.77 


4.46 


only 10% survives for T>100 K 


Cep E-mm 


Model 1 


1100 


1.48 


4.12 


0.71 






Model 2 


5000 


1.51 


0.86 


0.45 






Model 3 


1100 


0.91 


2.26 


1.19 


only 10% of CO evaporated 




Model 4 


1100 


1.42 


2.19 


1.12 


only 10% survives for T>100 K 


IC 1396 N 


Model 3 


1100 


2.95 


4.45 


>4.3 


only 10% of CO evaporated 




Model 4 


1100 


1.79 


4.35 


>4.0 


only 10% survives for T>100 K 


CBS 


Model 3 


1100 


1.63 


0.76 


0.31 


only 10% of CO evaporated 




Model 4 


1100 


1.24 


0.86 


0.51 


only 10% survives for T>100 K 


NGC 7129-FIRS 2 


Model 3 


1100 


0.97 


1.29 


2.05 


only 10% of CO evaporated 




Model 4 


1100 


1.39 


1.48 


1.73 


only 10% survives for T>100 K 



^ X is defined as x-'\l[^(^ model - hbsYM^- Note that this is an absolute error and the comparison among values in diff'erent sources is not 
straightforward. See Figs. 11 to 18. 



IC1396 error surface (RMS (K km/s)) IC1396 error surface (RMS (K km/s)) 



CB3 error surface (RMS (K km/s)) CB3 error surface (RMS (K km/s)) 




4000 6000 
Ro (AU) 




Ro (AU) 



Ro (AU) 



Fig. 9. The same as Fig. 8 for IC 1396 N. The solution with the 
lowest RMS, 0.36, is indicated by a cross. We consider accept- 
able solutions those with RMS less than 0.5. Contour levels are 
0.5, 1,2, and 5 K km s-^ 



et al. 2002), no atomic oxygen is included in the code because 
of the large uncertainties associated with its value (see e.g. Caux 
et al. 1999 and Melnick & Bergin 2005). This issue is discussed 
further at the end of this section. 

The abundances of the molecular ions (HCO^, N2H^, H3 , 
and all their deuterated forms) were calculated in terms of the 
instantaneous abundances of neutral species, assuming that the 
timescale for ion chemistry is much shorter than for freeze-out 
(Caselli et al. 2002). The rate coefficients are adopted from the 
UMIST database (http://www.udfa.net). For the proton-deuteron 
exchange reactions (such as H3 -h HD H2D^ -h H2), we used 
the rates measured by Gerlich et al. (2002), which better fit the 
deuterium fractionation in low-mass Class sources, as recently 
found by Emprechtinger et al. (2009). Hugo et al. (2009) have 
recently measured the proton-deuteron rate coefficients again. 



Fig. 10. The same as Fig. 8 for CB 3. The solution with the low- 
est RMS, 2.5, is indicated by a cross. Contour levels are 2.5, 4, 
6, and lOKkms'^ 



finding average values of total rates (for H3 and its deuterated 
isotopologues, the total rate refers to the average of multiple 
rates weighted according to the fraction of ortho and para forms; 
Sipila et al. 2010) about 4-5 times more than those derived by 
Gerlich et al. (2002, see also Sipila et al. 2010). Although this 
diflTerence could aflTect our results by increasing the deuterium 
fractionation by a similar factor (thus worsening the compari- 
son with observations), we decided not to include the new val- 
ues because the nuclear spin variants of all H3 isotopologues 
and H2 have not been distinguished in the current model. In 
fact, as Pagani et al. (1992) and Flower et al. (2004) showed, 
small temperature variations significantly alter the ortho-to-para 
ratio of H2, which in turn strongly aflTects the D-fractionation 
(ortho-H2 being more efficient than para-H2 in driving back 
the proton-deuteron exchange reactions thus decreasing the D- 
fractionation), even in the temperature regime between 9 and 
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20 K, where CO is mostly frozen onto dust grains (see also the 
discussion about the drop in ortho-H2D^ column density at tem- 
peratures above 10 K in Caselli et al. 2008). Such temperature 
variations are definitely present in the envelopes of intermediate- 
mass Class protostars (see also Emprechtinger et al. 2009), so 
that a simple increase in the rate coefficients without account- 
ing for nuclear spin variants and, in particular, the possible in- 
crease in the H2 ortho-to-para ratio may overestimate the D- 
fractionation calculated by our models. A more detailed chem- 
ical network is currently under development. The electron frac- 
tion has been computed using a simplified version of the reac- 
tion scheme of Umebayashi & Nakano (1990), where a generic 
molecular ion mH"^ is formed via proton transfer with H3 , and it 
is destroyed by dissociative recombination with electrons and re- 
combination on grain surfaces (using rates from Draine & Sutin 
1987). Dust grains follow a Mathis et al. (1977; MRN) size dis- 
tribution, but the minimum size has been increased to 5x10"^ 
cm to simulate possible dust coagulation, as in the best-fit model 
of the prestellar core LI 544 shown by Vastel et al. 2006 (see also 
Flower et al. 2005 and Bergin et al. 2006). The initial abundance 
of metals (assumed to freeze out with a rate similar to that of 
CO) is 10"^ (see McKee 1989). The cosmic ray ionization rate 
is fixed at ^=3x10"^^ s"^ (van der Tak & van Dishoeck 2000). 
The adopted CO binding energy is 1 100 K, a weighted mean of 
the CO binding energy on icy mantles (1180 K; CoUings et al. 
2003) and CO mantles (885 K; Oberg et al. 2005), assuming that 
solid water is about four times more abundant than solid CO. The 
models are computed for envelope lifetimes of 10^ yr, although 
solutions do not appreciably change after ^10^ yr. 

In Figs. 1 1 to 14, we show the fits for Cep E-mm. We adopt 
this protostar as a fiducial example because it is the one with the 
best spatial sampling of the envelope, its geometry looks spher- 
ical, and the contribution of the surrounding molecular cloud is 
negligible. It is impossible to reproduce the observations towards 
Cep E-mm using the standard chemical model. With the standard 
model, i.e., assuming a CO binding energy of -1100 K, we are 
able to reproduce the qualitative behavior of the C^^O and N2H^ 
emission (see Fig. 11), but the model reproduces the line inte- 
grated intensities poorly. The line intensity predictions are a fac- 
tor of 2-4 higher than the observations (see Fig. 11). Increasing 
the CO binding energy (e.g. assuming that a large fraction of CO 
is trapped in water ice) does not improve the fit (see model 2 in 
Table 6 and Fig. 12). 

The abundance profiles given in Table 5 suggest that the fit 
to the C^^O emission would improve by lowering the CO abun- 
dance in the inner part of the core. As pointed out there, this low 
CO abundance has some physical justification. One possibility is 
that a significant fraction of CO is converted into CH3OH on the 
grain surfaces before evaporation, in agreement with observa- 
tions of solid methanol along the line of sight towards embedded 
young stellar objects (e.g. Boogert et al. 2008). We mimic this 
situation with model 3. In this model, only 10% of the CO is re- 
leased back to the gas phase at the CO evaporation temperature. 
Another possibility is that CO is either destroyed by the stellar 
UV radiation and/or X-rays close to the star or transformed into 
more complex molecules via hot core chemistry. This case corre- 
sponds to our model 4, where only 10% of the CO survives when 
the temperature is higher than 100 K. The two models, 3 and 4, 
fit the C^^O and N2H^ emission better than the standard model 
(see Figs. 13 and 14), with model 3 being a slightly better fit to 
the Cep E-mm observations. Thus, we applied these most suc- 
cessful models, model 3 and 4, to all the other sources and have 
obtained reasonable fits (line integrated intensities fitted within 
a factor of 2) to the C^^O and N2H^ emission in all of them. For 



both IC 1396 N and Serpens-FIRS 1, model 4 gives a better fit, 
while model 3 is better for the rest. 

While we obtain reasonable fits for C^^O and N2H"^, our 
models do not succeed in reproducing the N2D^ data towards 
Serpens-FIRS 1 and IC 1396 N. In fact, the chemical models 
predict intensities higher by a factor of -10 than the observed 
intensities for these sources. All the considered models predict 
that the spatial distribution of N2D"^ is similar to that of N2H"^ 
and that the deuteration fraction [N2D+]/[N2H+]~0.01 in most 
of the envelope. This is not true for our low-deuterated sources, 
Serpens-FIRS 1 and IC 1396 N, in which the [N2D+]/[N2H+] 
ratio is a few 0.001. This discrepancy could have diflferent ori- 
gins. First of all, the models assume a spherical geometry. It is 
clear that bipolar outflows have excavated large cavities in these 
protostellar envelopes (see e.g. Fuente et al. 2009). The walls of 
these cavities are warmed by the stellar UV radiation and shocks. 
Moderate temperatures, UV radiation, and shocks would lower 
the abundance of N2H+ and change the [N2D+]/[N2H+] ratio. 
In fact, the large linewidths observed (~1.5 km s"^) are hard to 
maintain without shocks and dissipation. The role of UV radia- 
tion and shocks is expected to be greater in IM Class objects 
than in low-mass ones. 

This simple model ignores many other important parame- 
ters that could decrease the deuterium fractionation in warmer 
sources. These include (i) an increased abundance of atomic 
O in the gas phase, possibly coming from the release of water 
from the icy dust mantles. Atomic oxygen is an efficient de- 
struction partner for all the H3 isotopologues, and lowers the 
D-fractionation, as discussed by Caselli et al. (2002). (ii) The 
ortho-to-para H2 ratio, which in systems out of equilibrium (such 
as free-falling) can exceed the steady- state value by more than 
an order of magnitude (Flower et al. 2006, Pagani et al. 2009). 
A larger fraction of ortho-H2 leads to a lower D-fractionation, 
given that the more energetic ortho-H2 can more easily drive the 
proton-deuteron exchange reactions (e.g. H3 HD H2D"^ -\- 
H2) backward and reduce the H2D^/H3 abundance ratio (Gerlich 
et al. 2002). (iii) A higher ionization rate, f , which may be due 
to the presence of X-rays. A larger ^ implies a larger electron 
fraction and thus a higher dissociative recombination rates for 
molecular ions, including H3 isotopologues (e.g. Caselli et al. 
2008). (iv) The presence of small grains (in particular RAHs; see 
discussion in Caselli et al. 2008, Sect. 5.1) may arise from the 
interaction of outflow lobes and UV radiation with the molec- 
ular envelope. The associated shocks will partially destroy dust 
grains along the way (e.g. Jones et al. 1996; Caselli et al. 1997; 
Guillet et al. 2009) and increase the surface area for recombina- 
tion of molecular ions onto dust grains. 

Because of the low angular resolution of the present obser- 
vations, it is hard to disentangle the influences of the above pa- 
rameters. Both more detailed observations and more comprehen- 
sive models are needed for a deeper understanding of the chem- 
ical and physical evolution of the envelopes surrounding young 
intermediate-mass stars. 

8. Summary and conclusions 

We carried out a study of the CO depletion and N2H"^ deuteration 
in a sample of representative IM Class protostars. Our results 
can be summarized as follows. 

- We observed the millimeter lines of C^^O, C^^O, N2H+, and 
N2D"^ using the IRAM 30m telescope in a sample of 7 Class 
and 2 Class I IM stars. We have found a clear evolutionary 
trend that differentiates Class from Class I sources. While 



12 



T. Alonso-Albi et al.: Chemical study of intermediate-mass (IM) Class protostars. 

CepE 




Fig. 11. Left: results of our chemical model for Cep E-mm assuming a binding energy of 1100 K (standard value, model 1). The 
C^^O abundance is shown in black, N2H^ in red, and N2D^ in blue. Center: Comparison between the predicted C^^O 1^0 line 
intensities (continuous line) and the observed values (filled squares). Right: Same for N2H^ 1^0. 
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Fig. 12. The same as Fig. 1 1 for model 2. 
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Fig. 13. The same as Fig. 1 1 for model 3. 



the emission of the N2H^ 1^0 peaks towards the star posi- 
tion in Class protostars, it surrounds the FIR sources in the 
case of Class I stars. This occurs because the recently formed 
star has heated and disrupted the parent core in the case of 
Class I objects. The deuterium fractionation, R2, is low, be- 
low a few 0.001 in all Class I sources. There is, however, a 
wide dispersal in the values of R2 in Class sources ranging 
from a few 0.001 to a few 0.01. This at least two orders of 
magnitude greater than the elemental value in the interstellar 
medium, although a factor of 10-100 lower than in prestellar 
clumps. 

It is impossible, however, to establish an evolutionary trend 
among Class sources based on simple parameters such as 



the average CO depletion and the average N2H^ deuterium 
fractionation. This stems from the complexity of these re- 
gions (multiplicity) and the limited angular resolution of our 
observations, which prevents us from tracing the inner re- 
gions of the envelope. Interferometric observations are re- 
quired to provide a more precise picture of the evolutionary 
stage of these objects. 
- We used a radiative transfer code to derive the C^^O, N2H^, 
and N2D^ radial abundance profiles in 5 IM Class stars. In 
particular, we fit the C^^O 1^0 maps by assuming that the 
C^^O abundance decreases inwards within the protostellar 
envelope until the gas and dust reach the CO evaporation 
temperature, T^y. Our observational data are better fit with 
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Fig. 14. The same as Fig. 1 1 for model 4. 
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Fig. 15. The same as Fig. 1 1 for IC 1396 N and model 4. 
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Fig. 16. The same as Fig. 1 1 for CB 3 and model 3. 
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Fig. 17. The same as Fig. 1 1 for NGC 7129-FIRS 2 and model 3. 
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Serpens — Firs 1 




Fig. 18. The same as Fig. 1 1 for Serpens-FIRS 1 and model 3. 



values of Tgv~20-25 K, consistent with the binding energy 
of 1 100 K, which is the measured in the laboratory for a CO- 
CO matrix. 

- We determined the chemistry of the protostellar envelopes 
using the model by Caselli et al. (2002). A spherical envelope 
and steady- state chemical model cannot account for the ob- 
servations. We had to introduce modifications to better fit the 
C^^O and N2H^ maps. In particular the CO abundance in the 
inner envelope seems to be lower than the canonical value. 
This could be due to the conversion of CO into CH3OH on 
the grain surfaces, the photodissociation of CO by the stel- 
lar UV radiation, or even geometrical eff'ects. Likewise, we 
have problems fitting the low values of the deuterium frac- 
tionation (~ a few 0.001) measured for some Class IMs. 
Several explanations have been proposed to account for this 
discrepancy. 
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Appendix A: 

Below we discuss the details of the modeling for each individual 
source. 

A.1. Serpens-FIRS 1 

Serpens-FIRS 1, located near the center of the Serpens main 
core, is the most luminous object embedded in the cloud. Several 
continuum studies lead to its classification as a Class source 
with a bolometric luminosity estimated to range from 46 L© to 
84 Lo (Harvey et al. 1984; Casali et al. 1993; Hurt & Barsony 
1996; Larsson et al. 2000). The latest estimate of its luminos- 
ity (see Table 4) suggests that Serpens-FIRS 1 is on the low 
mass/IM borderline. Serpens-FIRS 1 drives a molecular out- 
flow that is oriented at a position angle of 50° (Rodriguez et al. 
1989). CIO modeled this source as a sphere with an outer radius 
of 5900 AU, 25.5" at the distance of Serpens. Convolving this 
size with the observational beams, one would expect a source 
diameter of 55", 57", and 53" for the C^^O 1^0, N2H+ 1^0, 
and N2D"^ 2^1 maps, respectively. In our maps, the emission 
is much more extended (>80"), showing that the dense molecu- 
lar cloud significantly contributes to the molecular emission (see 
also Fig. 6). To mimic the molecular cloud component, we in- 
creased the outer radius in the Crimier et al. profile. We assumed 
that the radius of the protostellar envelope is 15000 AU, and 
the dust density and temperature smoothly vary between the last 
point of the Crimier et al. profile (r=5900 AU, n=5xlO^ cm"-^, 
Tj=13 K) and the values assumed at r= 15000 AU, which are 
n^lxlO"^ cm"^ and Tj=10 K. This profile was used for our fit- 
ting. 

The integrated intensity map of C^^O cannot be fit with 
a constant abundance profile. Thus, we decided to assume a 
step function for the C^^O abundance: (i) a constant abun- 
dance, Xo, that is expected to be close to the canonical value 
(Xo=2±2xlO"^) for radii lower than a given radius, Rq and (ii) 
an abundance, Xi, that is expected to be <Xo for larger radii. 
Thus defined Rq, is the evaporation radius of CO. The values of 
Xo, Ro, and X\ were fit by the model. This approximation was 
still not sufficient to produce a good fit. We had to add another 
step, and two new variables, Ri and X2, in the C^^O abundance 
profile. The best fit is shown in Table 5. We have (i) a warm 
region (R<2000 AU) with an C^^O abundance -1.2x10"^ (ii) 
an intermediate layer with a high value for the C^^O depletion, 
f/)~20, and (iii) an external layer (R>6000 AU), which is essen- 
tially the molecular cloud component, in which the C^^O abun- 
dance is close to the canonical value again. 

We followed the same procedure for N2H^. In this case the 
emission extends to the NE and greatly differs from spherical 
symmetry. For this reason we masked the NE quadrant in our 
fitting (see Fig. 4). Again, we conclude that the N2H^ abun- 
dance has a standard value of a few 10"^^ in the inner region 
(R< 3000 AU) and in the molecular cloud, but decreases by at 
least a factor of 10 in the region between them. 

Serpens-FIRS 1 is one of the three sources where we 
were able to fit the N2D^ 2^1 emission. It is clear from 
the map morphology that the N2D"^ emission is dominated 
by the molecular cloud component. Similar to the case of 
N2H^, we masked the NE quadrant in our fitting. We found 
X(N2D^)= 1.9x10"^^ in the molecular cloud, and obtained an 
upper limit of X(N2D+)<2.0xlO"^^ for the protostellar core. 
Thus we have a deuterium fractionation of ~0.01 in the molec- 
ular cloud, and at least an order of magnitude lower in the core 
(see Fig. 7). 



A.2. Cep E 

Cep E-mm was cataloged as a Class protostar by Lefloch et 
al. (1996). Cep E-mm was observed with IRAM 30m (Lefloch 
et al. 1996; Chini et al. 2001), SCUBA (Chini et al. 2001), ISO 
(Froebrich et al. 2003), and Spitzer (Noriega-Crespo et al. 2005). 
All these studies confirm the Class status of Cep E-mm and 
constrain the source total mass and bolometric luminosity in 
the range of 7-25 M© and 80-120 L©, respectively. A bipolar 
molecular outflow, first reported by Fukui et al. (1989), is asso- 
ciated with Cep E-mm. The H2 and [Fell] study by Eisloffel et 
al. (1996) shows a quadrupolar outflow morphology suggesting 
that the driving source is a binary. 

CIO modeled this source as a sphere with an outer radius 
of 35800 AU, 49" at the distance of Cepheus (see Table 4). 
Convolving this size with the observational beams, one would 
expect a source diameter of ~108", ~110", and 100" for the 
C^^O 1^0, N2H+ 1^0, and N2D+ 2^1 maps, respectively. 
These sizes agree with those found in our maps. Of our sources, 
Cep E has the largest envelope diameter versus HPBW ratio and 
thus the best spatial sampling of the varying chemical conditions 
in its protostellar envelope (see Table 5). 

In this source, we fit the emission with a constant and close- 
to-standard abundance of 6x10"^ for radii below 3500 AU, and 
a power-law variation of the C^^O abundance for larger radii 
(see Table 5). The high value of fo, -10, would occur close 
to R=3500 AU. The same kind of profile was fit for the N2H+ 
abundance. In this case, the quality of the N2D"^ 2^1 line was 
not good enough to fit the abundance profile. We derived an av- 
eraged abundance across the envelope of ~7xl0"^^ by fitting the 
spectrum observed towards the center position. 

A3. IC 1396 N 

IC 1396 N is the globule associated with IRAS 21391+5802. 
It has strong submillimeter and millimeter continuum emission 
(Wilking et al. 1993; Sugitani et al. 2000; Codella et al. 2001), 
high density gas (Serabyn et al. 1993; Cesaroni et al. 1999; 
Codella et al. 2001; Beltran et al. 2004), and water maser emis- 
sion (FelH et al. 1992; Tofani et al. 1995; Patel et al. 2000). 
IC 1396 N is thus an active site of star formation. Using BIMA 
interferometric millimeter observations, Beltran et al. (2002) de- 
tected three sources (BIMA 1, BIMA 2, and BIMA 3) deeply 
embedded within the globule. Among the three, BIMA 2 has the 
strongest millimeter emission and is associated with an energetic 
E-W bipolar outflow (Codella et al. 2001 Beltran 2002, 2004). 
Recent studies by Neri et al. (2007) and Fuente et al. (2007) us- 
ing the Plateau de Bure interferometer show that BIMA 2 is itself 
a protocluster composed of 3 cores. 

CIO modeled BIMA 2 as a sphere with an outer radius of 
29600 AU, 39" at the distance of IC 1396 N (see Table 4). 
Together with Cep E-mm, it is our best-sampled protostellar en- 
velope (large envelope diameter to HPBW ratio). Interferometric 
observations published by Fuente et al. (2009) reveal that the 
outflow has already eroded a large biconical cavity in the molec- 
ular cloud. This is very likely the reason for the scarce C^^O 
1^0 and N2H"^ 1^0 emission to the east of the source. We have 
therefore masked this region to have a more accurate descrip- 
tion of the toroidal envelope (see Fig. 4). This masking does not 
affect our results significantly (more than a factor of 2 in the 
abundances). 

Similarly to Cep-E, the C^^O emission was better fit with a 
constant and close-to-standard abundance of 6x10"^ in the in- 
ner region (R<5500 AU) and a power-law variation of the C^^O 
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abundance for larger radii (see Table 5). The highest value of fo 
is ~5, and it would occur close to R=5500 AU. The N2H^ emis- 
sion was better fit with a two-step function, Xa^2^+ =4.2x10"^^ 
for R<10000 AU, 1.8x10-^^ for R>15000 AU, and LOxlO-^^ 
in between. The last step could be due to the vicinity of the 
sources BIMA 3 and BIMA 2 that heat the outer part of the en- 
velope. These sources are not considered in the (n-T) fit by CIO. 
In IC 1396 N, we have not detected the N2D+ 2^1 line. We de- 
rived an upper limit to the N2D^ abundance of <7xl0"^^. The 
lower value of fo and the non-detection of N2D^ is consistent 
with this source being a warmer and more evolved object. 

A.4. CB3 

The CB3 Bok globule is located at -2.5 kpc (Launhardt & 
Henning 1997; Wang et al. 1995). CB3-mm is the brightest 
millimeter source of the globule, first detected by Launhardt 
& Henning (1997) and subsequently observed in the sub- 
millimeter by Huard et al. (2000). Yun & Clemens (1994) de- 
tected a molecular bipolar outflow in CO, elongated in the 
NE-SW direction, associated with H2O masers (de Gregorio- 
Monsalvo et al. 2006). This outflow has been mapped in various 
molecular lines by Codella & Bachiller (1999), who concluded 
that it originates from CB3-mm. The same authors concluded 
that CB3-mm is probably a Class source. 

CB3 is diff'erent from all the other sources in our sample. In 
this object, the 24 fim sources are not spatially coincident with 
the 850 yum emission peak. CB3-mm hosts two 24 jim sources 
separated by ~12" and neither of them is spatially coincident 
with the column density peak. The column density peak, better 
traced by the 850 jim continuum emission, is located in between 
and almost equidistant from the two 24 //m sources (see Fig. 4). 
CIO modeled the SED of this source by adding up the flux of 
the two 24 yum sources. They fit the spatial distribution of the 
450 Jim and 850 jum maps as a sphere with an outer radius of 
103000 AU (i.e. 41") (see Table 4). The radial density power 
law in this source is steeper, ~2, than in the others, with low 
densities (~a few 10^ cm"^) in the outer envelope, > 50000 AU. 
This means that, although the protostellar envelope is very large, 
most of the emission comes from the inner 50000 AU. 

We fit the C^^O emission with a step function. The abun- 
dance is constant at 1.3x10"^ for radii less than 25000 AU, 
decreases to <3xlO"^ (//)>3), and then increases to 9.0x10"^ 
for R>60000 AU (see Table 5). The large C^^O abundance 
at large radii is very likely caused by the modeled extended 
low-density envelope providing a poor approximation. More 
likely, there are dense clumps immersed in a lower density 
cloud. The N2H^ emission is fit with a constant abundance 
of Xn^h+ =6. 5xlO~^^. The N2D^ emission is fit with a con- 
stant abundance of Xn^d+=^-OxIO~^^, implying an average 
[N2D+]/[N2H+]=0.05. 



tinuum emission from 25 to 2000 jam. Edwards & Snell (1983) 
detected a bipolar CO outflow associated with FIRS 2. The inter- 
ferometric study by Fuente et al. (2001) reveals that this outflow 
presents a quadrupolar morphology. In fact, the outflow seems 
to be the superposition of two flows, FIRS 2-out 1 and FIRS 2- 
out 2, likely associated with FIRS 2 and a more evolved star 
(FIRS 2-IR), respectively. Fuente et al. (2005a,b) carried out a 
complete chemical study of FIRS 2 providing the first detection 
of hot core in an IM Class 0. Based on all these studies, FIRS 2 
is considered the youngest IM object known at present. 

CIO modeled NGC 7129-FIRS 2 as a sphere with an outer 
radius of 18600 AU, 15" at the distance of NGC 7129 (see Table 
4). Since the source is very compact and barely resolved in the 
maps of the C^^O 1^0 and N2H^ 1^0 lines, we fit only an 
average abundance in for these species. The average CO deple- 
tion factor is fz)=3, consistent with the depletion factor found 
by Fuente et al. (2005a) based on H^^CO^ observations. The 
protostellar envelope is, however, resolved at the frequency of 
the N2D^ 3^2 line. Since previous interferometric observations 
published by Fuente et al. (2005a) reveal that the N2D^ emission 
is absent towards the hot core, we fitted the N2D"^ 3^2 emission 
assuming that it arises in a ring. The best fit is for a ring with in- 
ner radius of 9000 AU and an outer radius of 1 1000 AU. Within 
the ring the N2D^ abundance is ~4.5xl0"^\ and the deuterium 
fractionation, ~0.1, is the typical value found in prestellar cores. 



A.5. NGC 7129-FIRS 2 

NGC 7129 is a reflection nebula located in a complex and active 
star-forming site at a distance of 1250 pc (Hartigan & Lada 1985; 
Miranda et al. 1993). NGC 7129 FIRS 2 is not detected at opti- 
cal or near-infrared wavelengths. Its position is spatially coinci- 
dent with a ^^CO column density peak (Bechis et al. 1978) and 
a high-density NH3 cloudlet (Guesten & Marcaide 1986), and it 
is close to an H2O maser (Rodriguez et al. 1980). NGC7129- 
FIRS 2 has been classified as a Class IM protostar by Eiroa et 
al. (1998), who carried out a multi- wavelength study of the con- 



